Phenotypic responses to temperature in the ciliate Tetrahymena thermophila

Study summary

In this study, we grew four lineages of the ciliate Tetrahymena thermophila at 38 °C, a highly stressful temperature. We have monitored population abundances and two important morphological traits, cell size and cell shape, throughout the experiment. In the middle of the experiment, we moved half of the populations back to the control temperature of 20 °C, to evaluate possible costs related to the exposure to this high temperature environment.

Temperature range pilot experiment

We performed a preliminary experiment to determine the temperature range in which this strain of T. thermophila (1630/1U) can grow.

Load the growth data across 10 temperatures.

library(tidyverse)
library(ggpubr)
library(here)
library(patchwork)

temperature_range_data <- read_csv(here("1-data/TempChoice_abundance_dryad.csv"))

Create a plot of population abundance through time for each population replicate, coloring them by growth temperature.

Figure 1 Population dynamics of T. thermophila strain 1630/1U growing in ten different temperatures.

Each line represents one replicate population and the colors indicate the temperature in which the population was grown.

Temperature_ranges_plot <- ggplot(data = temperature_range_data, aes(x = day, y = mean_density, group = population, color = as.factor(temperature))) +
  geom_point() +
  geom_line() +
  scale_y_log10() +
  xlab("Day") +
  ylab("Population density (cells/ml)") +
  scale_color_manual(name = "Temperature (°C)", 
                     values = c('#053061', '#2166ac', '#4393c3', '#92c5de',
                                '#d1e5f0','#fddbc7','#f4a582','#d6604d','#b2182b', '#67001f')) +
  theme(text = element_text(size = 15),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.position = "bottom",
        plot.margin = unit(c(1,1,1.5,1.2),"cm"))

Temperature_ranges_plot

Temperature experiment

Figure 2B Experimental design of the temperature experiment. In the schematic of the experiment, each bottle represents one replicate batch culture, and the colors indicate the temperature in which the culture was grown.

Load the population abundance and the morphology data

# load the abundance data
abundance_data <- read_csv(here("1-data/TempAdapt_abundances_dryad.csv"))

# load the morphology data
morphology_data <- read_csv(here("1-data/TempAdapt_morphology_dryad.csv"))

Calculate the total number of measured cells and the mean number of measured cells per population per day

# total number of cells
unique(morphology_data$id) %>% length()

## [1] 94344

# mean number of cells plus standard error
cells_summary <- morphology_data %>% group_by(population, date) %>% summarise(measured_cells = n())

cells_summary %>% ungroup() %>% summarise(mean_measured_cells = mean(measured_cells), se_measured_cells = sd(measured_cells)/sqrt(n()))

## # A tibble: 1 x 2
##   mean_measured_cells se_measured_cells
##                 <dbl>             <dbl>
## 1                513.              26.9

Calculate mean, variance, standard deviation (sd) and standard error (se) of cell size and cell shape per population per day

# calculate mean, variance and sd of area and AR for each population for each day

morphology_data_ts <- morphology_data %>% group_by(population, lineage, batch, date, day, temperature) %>%
  summarise(mean_area_day = mean(mean_area), sd_area_day = sd(mean_area), 
            se_area_day = sd(mean_area)/sqrt(n()), var_area_day = var(mean_area),
            mean_AR_day = mean(mean_ar), sd_AR_day = sd(mean_ar), 
            se_AR_day = sd(mean_ar)/sqrt(n()))

Calculate coefficient of variation of cell area

# calculate coefficient of variation of cell area
morphology_data_ts <- morphology_data_ts %>% mutate(cv_area_day = sd_area_day/mean_area_day)

Merge morphology and abundances dataset

# merge both datasets
morphology_abundance_data <- left_join(abundance_data, morphology_data_ts)

# add name column for plotting (lineage + temperature)
morphology_abundance_data <- morphology_abundance_data %>% mutate(name = paste0(lineage, "_", temperature))

Calculate number of generations in each batch

# number of generations per batch ####
generations <- morphology_abundance_data %>% group_by(population, lineage, batch,replicate, temperature) %>% 
  summarise(maxi_density = max(mean_density), min_density = min(mean_density),
            number_generations = (log(max(mean_density)) - log(min(mean_density)))/log(2))

generations_summary <- generations %>% group_by(batch, temperature) %>% summarise(mean_generations_batch = mean(number_generations))

# create a tibble with the number of generations per batch at each temperature
generations_summary <- generations_summary %>% group_by(batch) %>% pivot_wider(names_from =  temperature, values_from = mean_generations_batch, names_prefix = "generations_")

# round the number of generations
generations_summary <- generations_summary %>% mutate(text_plot = case_when(
  batch == 1 ~ paste0(round(generations_20, digits = 0), " gen."),
  batch == 2 | batch == 3 ~ paste0(round(generations_38, digits = 0), " gen."),
  batch > 3 ~ paste0(round(generations_20, digits = 0), " gen. 20°C \n", round(generations_38, digits = 0), " gen. 38°C")))

Plot time series of population abundance

# define size of small and large fonts for the plot
small_font <- 10
large_font <- 13
  
plot_abundances <- ggplot(data = morphology_abundance_data) +
  geom_point(aes(x = batch_day, y = mean_density, group = population, colour = name), size = 1) + 
  geom_line(aes(x = batch_day, y = mean_density, group = population, colour = name)) +
  geom_label(data = generations_summary, aes(x = 9, y = 500, label = text_plot), size = 2.1) +
  scale_y_log10() +
  ylab("Population density\n(cells/ml)") +
  xlab("") +
  scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "red", "blue"), 
                     guide = F, name = "Population and\n temperature", 
                     labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                            "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                            "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C", 
                                            "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  scale_x_continuous(breaks = c(5, 10)) +
  theme(text = element_text(size = large_font),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.position = "bottom",
        plot.margin = unit(c(0, 1.5, 0 , 1.5 ),"cm")) +
  facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2", 
                                                          "3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5"))) 

Plot time series of each morphology measurement

# minimum and maximum cell area in the control
area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>% 
  select(mean_area_day) %>% min()
area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>% 
  select(mean_area_day) %>% max()

# plot cell area
plot_cell_area <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = mean_area_day, 
                                                               group = population, col = name)) +
  geom_point(size = 1) +
  geom_hline(aes(yintercept = area_min), linetype = "dashed") +
  geom_hline(aes(yintercept = area_max), linetype = "dashed") +
  geom_errorbar(aes(ymin = mean_area_day - se_area_day, ymax = mean_area_day + se_area_day, width = 0.9), size = 0.2) +
  geom_line() +
  xlab("") +
  ylab("Cell size (um2)") +
  scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
                     name = "Population and\n temperature", 
                     labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                            "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                            "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                            "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  scale_x_continuous(breaks = c(5, 10)) +
  theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.position = "bottom", strip.text = element_blank(), 
        plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
  facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2", 
                                                         "3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))


# minimum and maximum cv of cell area in the control
cv_area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>% 
  select(cv_area_day) %>% min()
cv_area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>% 
  select(cv_area_day) %>% max()

# plot coefficient of variation of cell area
# multiply cv per 100 to obtain percentage
plot_cv_cell_area <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = cv_area_day*100, 
                                                                  group = population, col = name)) +
  geom_point(size = 1) +
  geom_hline(aes(yintercept = cv_area_min*100), linetype = "dashed") +
  geom_hline(aes(yintercept = cv_area_max*100), linetype = "dashed") +
  geom_line() +
  xlab("") +
  ylab("Coefficient of variation\n of cell size (%)") +
  scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
                     name = "Population and\n temperature", 
                     labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                            "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                            "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                            "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  scale_x_continuous(breaks = c(5, 10)) +
  scale_y_continuous(breaks = c(15, 30, 45, 60)) +
  theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.position = "bottom", strip.text = element_blank(), 
        plot.margin = unit(c(0, 1.5, 0, 1.5 ),"cm")) +
  facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2", 
                                                         "3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))


# minimum and maximum variance of cell area in the control
var_area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>% 
  select(var_area_day) %>% min()
var_area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>% 
  select(var_area_day) %>% max()

# plot variance of cell area
plot_var_cell_area <-  ggplot(data = morphology_abundance_data, aes(x = batch_day, y = var_area_day, group = population, col = name)) +
  geom_point(size = 1) +
  geom_hline(aes(yintercept = var_area_min), linetype = "dashed") +
  geom_hline(aes(yintercept = var_area_max), linetype = "dashed") +
  geom_line() +
  xlab("Day") +
  ylab("Variance of\n cell size")  +
  scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
                     name = "Population and\n temperature", 
                     labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                            "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                            "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                            "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  scale_x_continuous(breaks = c(5, 10)) +
  theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.position = "bottom", strip.text = element_blank(), 
        plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
  facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2", 
                                                         "3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum cell aspect ratio in the control
AR_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>% 
  select(mean_AR_day) %>% min()
AR_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>% 
  select(mean_AR_day) %>% max()

# plot cell aspect ratio
plot_cell_AR <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = mean_AR_day, 
                                                             group = population, col = name)) +
  geom_point(size = 1) +
  geom_hline(aes(yintercept = AR_min), linetype = "dashed") +
  geom_hline(aes(yintercept = AR_max), linetype = "dashed") +
  geom_errorbar(aes(ymin = mean_AR_day - se_AR_day, ymax = mean_AR_day + se_AR_day, width = 0.9), size = 0.2) +
  geom_line() +
  xlab("") +
  ylab("Cell shape")  +
  scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
                     name = "Population and\n temperature", 
                     labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                            "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                            "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                            "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  scale_x_continuous(breaks = c(5, 10)) +
  theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.position = "bottom", strip.text = element_blank(), 
        plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
  facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2", 
                                                         "3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))

Figure 3 Population dynamics and morphological traits of each T. thermophila population during the temperature experiment.

Population density (A), mean cell size (B), variance of cell size (C), coefficient of variation of cell size (D) and mean cell shape (E) are shown for each population and for each batch separately. Minimum number of generations that took place in each batch is shown in boxes in plot A. Error bars indicate standard errors of means for cell size and cell shape (B and E). The colors indicate the temperature in which the population was grown, and the shade represents the population replicate. Dashed lines mark the range of observed values at the control temperature (20 °C) in the first batch of the experiment.

# all plots in one figure
all_time_series <- ggarrange(plot_abundances, plot_cell_area, plot_var_cell_area, plot_cv_cell_area, plot_cell_AR,
                             nrow = 5, common.legend = T, legend = "bottom", align = "hv", 
                             labels = "AUTO")

# print the plot
all_time_series

Estimate growth rate and lag phase in each population using the Gompertz model.

# load the package growthrates
library(growthrates)

# select only relevant data for growth model (population, population density and temperature)
# exclude populations 32, 141 and 341, since model will be run separetely for them
data_subset1 <- morphology_abundance_data %>%
  group_by(population, lineage, batch, replicate, temperature) %>% 
  select(population, mean_density, batch_day) %>%
  filter(!population %in% c("32", "141", "341"))

# set parameters to start the model (p), lower boundary and upper boundary
# y0 =  initial population density, mumax = maximum growth rate, K = carrying capacity, lambda = lag phase

p     <- c(y0 = 100,   mumax = 1,   K = 1000000, lambda = 5)
lower <- c(y0 = 0,      mumax = 0,   K = 100000,  lambda = 0)
upper <- c(y0 = 5000,   mumax = 15,  K = 5000000, lambda = 8)

# run gompertz model with lag phase estimation for each population separately
subset1_gompertz3 <- all_growthmodels(
  mean_density ~ grow_gompertz3(batch_day, parms) | population,
  data = data_subset1,
  p = p, lower = lower, upper = upper,
  log = "y", ncores = 4)

Check model fit and coeficients.

# check the rsquared of the models

# r squared of the model
rsquared(subset1_gompertz3) > 0.9

##    11    12    21    22    31    41    42   131   132   142   151   152 
##  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
##   231   232   241   242   251   252   331   332   342   351   352   431 
##  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE 
##   432   441   442   451   452 
##  TRUE  TRUE  TRUE  TRUE  TRUE

rsquared(subset1_gompertz3) # all above 0.90, expect pop. 241, with 0.80

##        11        12        21        22        31        41        42 
## 0.9883432 0.9178892 0.9933680 0.9998310 0.9937436 0.9809305 0.9981437 
##       131       132       142       151       152       231       232 
## 0.9998194 0.9998437 0.9998797 0.9914330 0.9998212 0.9999818 0.9999183 
##       241       242       251       252       331       332       342 
## 0.8075650 0.9992203 0.9998438 0.9998244 0.9998084 0.9539876 0.9999996 
##       351       352       431       432       441       442       451 
## 0.9477612 0.9993824 0.9992420 0.9996416 0.9997832 0.9999890 0.9997705 
##       452 
## 0.9994084

# coefficients of the model
coef(subset1_gompertz3)

##             y0     mumax         K       lambda
## 11  4999.99925  1.197897 1129221.8 5.426150e-06
## 12   490.11662  3.939040  156097.9 7.999999e+00
## 21  2822.39743  1.342365 1261845.6 1.789159e-06
## 22   762.82078  1.764727  584705.5 7.973149e+00
## 31   627.35002  1.906846 1173617.9 2.705005e-06
## 41  4999.99917  1.173631 1373269.9 4.339097e-06
## 42  1736.92242  1.354624  100000.0 7.647302e+00
## 131 1654.64387  1.676791 1081410.9 7.102917e+00
## 132 1589.40401  1.251949 3011394.3 6.268673e+00
## 142 1452.04283  1.144168 4999513.5 4.818932e-01
## 151   19.72001  4.887823  248032.6 4.603856e-01
## 152 1335.03526  1.412242 2279297.1 2.774922e-01
## 231 1426.42616  2.897612  504515.5 7.408923e+00
## 232  924.84138  2.550602 4997991.7 7.924815e+00
## 241 3241.76987 14.999383  226100.6 2.869473e+00
## 242 4999.99967  2.225997  743961.4 2.876162e+00
## 251  399.02617  2.094700 4685750.6 1.448947e+00
## 252  855.05594  2.145062  924045.8 9.043278e-01
## 331 2093.40965  1.539841 4456116.9 3.773430e+00
## 332 4313.96423 14.999743  264803.7 4.813319e+00
## 342  806.37213  1.727175 1957003.6 6.825973e-01
## 351 2024.29552 14.760898  201141.4 1.765625e+00
## 352 1050.07698  1.509381 1850540.4 2.545231e-01
## 431  107.09779  2.026098  274606.3 3.168215e+00
## 432 1341.19534  2.204486  137174.9 4.250084e+00
## 441 1618.63891  1.823084  167147.9 2.417790e+00
## 442 1558.57363  2.097457  767992.8 1.765452e+00
## 451 2814.32870  3.254367  264303.9 1.972106e+00
## 452  583.12852  2.192975  808572.3 7.190942e-01

# plot the model, y in log scale
par(mfrow = c(8, 4))
par(mar = c(2.5, 4, 2, 1))
plot(subset1_gompertz3)
plot(subset1_gompertz3, log = "y")

Gompertz model for populations that need specific parameters

# populations 32, 141, 341 presented suboptimal fittings with the above parameters
# model these populations separately

# population 32, remove last day since population is already crashing
pop32 <- morphology_abundance_data %>%
  group_by(population, lineage, batch, replicate, temperature) %>% 
  select(population, mean_density, batch_day) %>% 
  filter(population == 32 & batch_day < 12)

## Adding missing grouping variables: `lineage`, `batch`, `replicate`, `temperature`

# same parameters as previous populations
pop32_gompertz <- fit_growthmodel(FUN = grow_gompertz3, p = p,
                                  pop32$batch_day, pop32$mean_density,
                                  lower = lower, upper = upper)

# # population 141
pop141 <- morphology_abundance_data %>%
  group_by(population, lineage, batch, replicate, temperature) %>%
  select(population, mean_density, batch_day) %>%
  filter(population == 141 & batch_day < 6)

## Adding missing grouping variables: `lineage`, `batch`, `replicate`, `temperature`

# use new parameters
p2     <- c(y0 = 500,    mumax = 1, K = 1000000, lambda = 0.1)
lower2 <- c(y0 = 0,      mumax = 0,   K = 100000,   lambda = 0)
upper2 <- c(y0 = 10000,   mumax = 30,  K = 5000000, lambda = 4)

pop141_gompertz <- fit_growthmodel(FUN = grow_gompertz3, p = p,
                                   pop141$batch_day, pop141$mean_density,
                                   lower = lower, upper = upper)


# population 341, remove final day that population is crashing
pop341 <- morphology_abundance_data %>%
  group_by(population, lineage, batch, replicate, temperature) %>% 
  select(population, mean_density, batch_day) %>%
  filter(population == 341 & batch_day <= 5)

## Adding missing grouping variables: `lineage`, `batch`, `replicate`, `temperature`

# rsquared of the model is very low, manually estimate growth rate and lag phase

pop341_day0 <- pop341 %>% filter(batch_day == 0)
pop341_day3 <- pop341 %>% filter(batch_day == 3)

# calculate maximum growth rate as the log difference in population density divided by the time difference
mumax <- (log10(pop341_day3$mean_density) - log10(pop341_day0$mean_density))/(pop341_day3$batch_day - pop341_day0$batch_day)
lambda <- 0 # define lag phase as 0
y0 <- NA
K <- NA
population <- 341

pop341_manual_growth <- tibble(mumax, lambda, y0, K, population)


remove(mumax, lambda, y0, K, population)

Check model fit and coeficients.

# pop 32
rsquared(pop32_gompertz)

## [1] 0.8966222

coef(pop32_gompertz)

##           y0        mumax            K       lambda 
## 2.969880e+03 7.143717e+00 1.681196e+05 7.999997e+00

par(mfrow = c(1,2))
plot(pop32_gompertz)
plot(pop32_gompertz, log = "y")

# pop 141

rsquared(pop141_gompertz)

## [1] 0.9996588

coef(pop141_gompertz)

##           y0        mumax            K       lambda 
## 2.047749e+03 3.552567e+00 2.361080e+05 2.279727e+00

par(mfrow = c(1,2))
plot(pop141_gompertz)
plot(pop141_gompertz, log = "y")

Obtain model predictions

# try to get the model predictions for all gompertz models

# new time variable with the maximum batch_day across all populations (12)
time <- data.frame(time = seq(from = 0, to = 12, by = 1))

# predict function for each problematic population, add name of the population to the tibble
pop32_growth_prediction <- pop32_gompertz %>% predict(newdata = time) %>% as_tibble() %>% dplyr::rename(batch_day = time, mean_density_predict = y) %>% mutate(population = "32")
pop141_growth_prediction <- pop141_gompertz %>% predict(newdata = time) %>% as_tibble() %>% dplyr::rename(batch_day = time, mean_density_predict = y) %>% mutate(population = "141")

# predict function for each the remaining populations
subset1_growth_prediction <- predict(subset1_gompertz3, newdata = time)
names(subset1_growth_prediction)

##  [1] "11"  "12"  "21"  "22"  "31"  "41"  "42"  "131" "132" "142" "151"
## [12] "152" "231" "232" "241" "242" "251" "252" "331" "332" "342" "351"
## [23] "352" "431" "432" "441" "442" "451" "452"

# transform list into data frame with population as a column
subset1_growth_prediction <- map_df(subset1_growth_prediction, ~as.data.frame(.x), .id = "population")

# rename the columns
subset1_growth_prediction <- subset1_growth_prediction %>% dplyr::rename(batch_day = time, mean_density_predict = y)

# merge all models
prediction_gompertz <- bind_rows(subset1_growth_prediction, pop32_growth_prediction, pop141_growth_prediction)
prediction_gompertz$population <- as.numeric(prediction_gompertz$population)

morphology_abundance_data_gompertz <- morphology_abundance_data %>% left_join(prediction_gompertz)

Merge all growth model estimates

# create a dataset with the coeficients of the model
coef_gompertz <- coef(subset1_gompertz3) %>% as_tibble()
coef_gompertz <- cbind(coef_gompertz, population = as.double(names(subset1_gompertz3)))

coef_gompertz <- coef_gompertz %>% rbind(c(coef(pop32_gompertz),  32))
coef_gompertz <- coef_gompertz %>% rbind(c(coef(pop141_gompertz),  141))
coef_gompertz <- coef_gompertz %>% rbind(pop341_manual_growth)

# select metadata from the original dataset
metadata <- morphology_abundance_data %>% select(population, lineage, batch, replicate, temperature) %>% distinct()

# add metadata to the growth model coeficients data
coef_gompertz <- coef_gompertz %>% left_join(metadata)

## Joining, by = "population"

# add name column for plotting (lineage + temperature)
coef_gompertz <- coef_gompertz %>% mutate(name = paste0(lineage, "_", temperature))

Figure 4 Growth dynamics of T. thermophila during the temperature experiment.

Gompertz model was used to estimate demographic parameters. (A) Points show population abundances and lines show the fitting of the model. No model fit is shown for population 3 in batch 4 at 38 °C since growth parameters were manually calculated (see methods). Maximum growth rate (B) and lag phase (C) of each T. thermophila population per batch culture. In all plots, the colors indicate the temperature in which the population was grown, and the shades represent the population replicate. The lag phases in batch 1 are estimated as smaller than 1x10-5 and are not visible in the plot.

# define font sizes
small_font <- 10
large_font <- 13

# plot the lag phase length (lambda)
plot_lambda <- ggplot(coef_gompertz, aes(x = batch, y = lambda, fill = name, group = population)) +
  geom_col(position = "dodge") +
  xlab("") +
  ylab("Lag phase (days)") +
  scale_fill_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
                    name = "Population and\n temperature",
                    labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                           "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                           "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                           "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  guides(fill = guide_legend(title.theme = element_text(size = small_font), label.theme = element_text(size = small_font))) +
  theme(text = element_text(size = large_font),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_blank(),
        legend.position = "right",
        #plot.margin = margin(t = 10, r = 10, b = 10, l = 10),
        axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
        axis.text.x = element_blank(), axis.ticks.x = element_blank()) + #,
  facet_grid(~batch, scales = "free_x", 
             labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2", 
                                      "3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))


# plot the maximum growth rate (mumax)
plot_mumax <- ggplot(coef_gompertz, aes(x = batch, y = mumax, fill = name, group = population)) +
  geom_col(position = "dodge") +
  xlab("") +
  ylab(bquote("Max. growth rate"~(day^-1))) +
  scale_fill_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
                    name = "Population and\n temperature", guide = F,
                    labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                           "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                           "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                           "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  #guides(fill = guide_legend(title.theme = element_text(size = small_font), label.theme = element_text(size = small_font))) +
  theme(text = element_text(size = large_font),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_blank(),
        legend.position = "bottom",
        #plot.margin = margin(t = 10, r = 10, b = 10, l = 10),
        axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
        axis.text.x = element_blank(), axis.ticks.x = element_blank()) + #,
  facet_grid(~batch, scales = "free_x", 
             labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2", 
                                      "3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))



gompertz_fit_plot <-  ggplot(data = morphology_abundance_data_gompertz) +
  geom_point(aes(x = batch_day, y = mean_density, group = population, colour = name), size = 1.5) +
  geom_line(aes(x = batch_day, y = mean_density_predict, group = population, colour = name)) +
  scale_y_log10() +
  ylab("Population density (cells/ml)\n") +
  xlab("Day\n") +
  scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "red", "blue"),
                     name = "Population and\n temperature", guide = F,
                     labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                            "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                            "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                            "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  scale_x_continuous(breaks = c(5, 10)) +
  theme(text = element_text(size = large_font),
        panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA),
        axis.line = element_line(colour = "black"),
        legend.position = "bottom",
        plot.margin = unit(c(0, 1.5, 0 , 1.5 ),"cm")) +
  facet_grid(paste0("Batch ", batch) ~ paste0("Population ", lineage))


# layout of the final figure
layout <- "
AA#
AA#
AA#
BBD
CCD"

# create final figure with three facets
plot_growth <- gompertz_fit_plot + plot_mumax + plot_lambda + guide_area() + plot_layout(guides = "collect", design = layout) + plot_annotation(tag_levels = 'A')

# print figure
plot_growth

Temperature response modelling

Calculate the change in growth rate, lag phase and in morphological traits compared to the control.

The control is the mean of the four populations in batch 1, growing at 20 °C.

# create variable merging batch, temperature and replicate; select only used variables
growth_response_gompertz <- coef_gompertz %>% ungroup() %>% 
  mutate(batch_temp = paste0("b", batch, "_", temperature, "_", replicate)) %>%
  select(-population, -batch, - replicate, - temperature, -name)  

# calculate % change in comparison to the same population lineage in batch 1
growth_response_gompertz <- growth_response_gompertz %>% select(-y0, -K) %>% pivot_longer(cols = 1:2, names_to = "data_type")
  
growth_response_gompertz <- growth_response_gompertz %>% group_by(lineage, data_type) %>% pivot_wider(names_from = batch_temp, values_from = value)

# calculate the difference in maximum growth rate and in lag phase duration in absolute values, not the % change since numbers are too large and interpretation is difficult
growth_response_gompertz <- growth_response_gompertz %>% 
  mutate(week2 = b2_38_NA - b1_20_NA,
         week3_1 = b3_38_1 - b1_20_NA,
         week3_2 = b3_38_2 - b1_20_NA,
         week4 = b4_38_1 - b1_20_NA,
         week5 = b5_38_1 - b1_20_NA,
         week4_cold = b4_20_2 - b1_20_NA,
         week5_cold = b5_20_2 - b1_20_NA)


# tidy the data
growth_response_gompertz <- growth_response_gompertz %>% select(1, 2, 11:17) %>% 
  gather(key = "comparison", value = "percent_batch1", 3:9)


# repeat the same calculations with the morphological traits

# calculate mean cell area and cell aspect ratio in each batch
morphology_data_batch <- morphology_abundance_data %>% drop_na(mean_area_day) %>%
  group_by(population, lineage, batch, replicate, temperature) %>%
  summarise(mean_area_batch = mean(mean_area_day), cv_area_batch = mean(cv_area_day), var_area_batch = mean(var_area_day),
            mean_AR_batch = mean(mean_AR_day))

# create variable combining batch, temperature and replicate; select only used variables
morphology_response <- morphology_data_batch %>% ungroup() %>% 
  mutate(batch_temp = paste0("b",batch, "_", temperature, "_", replicate)) %>%
  select(lineage, batch_temp, mean_area_batch, mean_AR_batch, cv_area_batch, var_area_batch) 

morphology_response <- morphology_response %>% gather(key = "data_type", value = "value", 3:6)

morphology_response <- morphology_response %>% group_by(lineage, data_type) %>% spread(key = batch_temp, value = value)

# calculate % change in comparison to the same population lineage in batch 1
morphology_response <- morphology_response %>% group_by(lineage, data_type) %>% 
  mutate(week2 = ((b2_38_NA - b1_20_NA) * 100) / b1_20_NA,
         week3_1 = ((b3_38_1 - b1_20_NA)  * 100) / b1_20_NA,
         week3_2 = ((b3_38_2 - b1_20_NA)  * 100) / b1_20_NA,
         week4 = ((b4_38_1 - b1_20_NA)  * 100) / b1_20_NA,
         week5 = ((b5_38_1 - b1_20_NA)  * 100) / b1_20_NA,
         week4_cold = ((b4_20_2 - b1_20_NA)  * 100) / b1_20_NA,
         week5_cold = ((b5_20_2 - b1_20_NA)  * 100) / b1_20_NA)

morphology_response <- morphology_response %>% select(1, 2, 11:17) %>% 
  gather(key = "comparison", value = "percent_batch1", 3:9)

# merge growth and lag phase with the morphology data
response_data <- rbind.data.frame(growth_response_gompertz, morphology_response)

# add batch information
response_data <- response_data %>% ungroup() %>% mutate(batch = case_when(
  comparison == "week2" ~ 2,
  comparison == "week3_1" ~ 3,
  comparison == "week3_2" ~ 3,
  comparison == "week4" ~ 4,
  comparison == "week5" ~ 5,
  comparison == "week4_cold" ~ 4,
  comparison == "week5_cold" ~ 5
))

# add temperature information
response_data <- response_data %>% mutate(temperature = case_when(
  comparison == "week2" ~ 38,
  comparison == "week3_1" ~ 38,
  comparison == "week3_2" ~ 38,
  comparison == "week4" ~ 38,
  comparison == "week5" ~ 38,
  comparison == "week4_cold" ~ 20,
  comparison == "week5_cold" ~ 20
))

# add variable for lineage and temperature color
response_data <- response_data %>% mutate(name = paste0(lineage, "_", temperature))


# change data type into factor and change levels for a nicer looking plot
response_data$data_type <- factor(response_data$data_type, levels = c("mumax", "lambda", 
                                                                              "mean_area_batch", "cv_area_batch", "var_area_batch", 
                                                                              "mean_AR_batch"))

# transform lineage into categorical variable
response_data$lineage <- as.character(response_data$lineage) 

# look at the data
ggplot(data = response_data, aes(x = batch, y = percent_batch1, color = name)) +
  geom_point(alpha = 0.8, size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "#878787"),
                     name = "Population and\n temperature", 
                     labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                            "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                            "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                            "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  ylab("Change compared to control (%)") +
  xlab("Batch culture") +
  guides(fill = guide_legend(title.theme = element_text(size = 15), 
                             label.theme = element_text(size = 15))) +
  theme(text = element_text(size = 15), 
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_blank(),
        legend.position = "bottom") +
  facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
                                                            "lambda" = "Lag phase",
                                                            "cv_area_batch" = "CV cell size",
                                                            "mean_AR_batch" = "Cell shape",
                                                            "mean_area_batch" = "Cell size",
                                                            "var_area_batch" = "Variance cell size")), scales = "free_y") 

Model the response with Bayesian mixed-effects models using the package MCMCglmm.

1) Response of the populations that grew at 38 °C.

# load the libraries
library(MCMCglmm)
library(coda)
library(purrr)

set.seed(1)

# select the populations that experienced 38 °C
response_data_38 <- response_data %>% filter(temperature == 38)

# create variable for time difference since warming started
response_data_38 <- response_data_38 %>% mutate(batches_38 = batch - 2)

# create a list with the prior information, using an uninformative prior
prior_1 <- list(R = list(V = 1, nu = 0.002), 
                G = list(G1 = list(V = 1, nu = 0.002)))

# run one separate model for each of the 5 response variables
models_38 <- response_data_38 %>% group_by(data_type) %>% nest() %>%
  mutate(mixed_model = map(data, ~MCMCglmm(fixed = percent_batch1 ~ batches_38, 
                                           random = ~lineage, 
                                           data = ., 
                                           family = "gaussian",
                                           prior = prior_1, 
                                           nitt = 2000000, 
                                           burnin = 30000, 
                                           thin = 1000, verbose = F)))


# create new data
new_data <- data.frame(expand.grid(lineage = c("1"), batches_38 = c(0:3), 
                       percent_batch1 = 0))

# get confidence interval from each model
models_38 <- models_38 %>% mutate(predictions = map(mixed_model, ~predict.MCMCglmm(object = .,
                                                               newdata = new_data,
                                                               interval = "confidence")))
models_38 <- models_38 %>% mutate(predictions = map(.x = predictions, ~as_tibble(.)))

# merge new data and model fit for plotting
models_38 <- models_38 %>% mutate(predictions = map(.x = predictions, ~cbind(new_data, .))) 

# store models in a separate object
model_38_predictions <- unnest(models_38, predictions)

# plot the data and the model 
response_38_plot <- ggplot(data = response_data_38, aes(x = batch, y = percent_batch1, color = name)) +
  geom_point(aes(color = as.factor(lineage)), alpha = 0.8, size = 2) +
  geom_smooth(data = model_38_predictions, aes(x = batches_38 + 2, y = fit, ymin = lwr, ymax = upr), colour = "black", stat = "identity") +
  geom_hline(yintercept  = 0, linetype = "dashed") +
  #scale_y_continuous(limits = c(-95, 95)) +
  xlab("Batch") +
  ylab("Change in comparison to control (%)") +
  scale_color_manual(values = c("#efd0d4", "#d88b95","#c14655","#b2182b"),
                     name = "Population", labels = as_labeller(c( "1" = "Pop. 1",
                                                                  "2" = "Pop. 2",
                                                                  "3" = "Pop. 3",
                                                                  "4" = "Pop. 4"))) +
  theme(text = element_text(size = 15), plot.title = element_text(hjust = 0.5, margin = unit(c(0.5,0.5,0.5,0.5), "lines")),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_blank(),
        legend.position = "bottom") +
  facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
                                                            "lambda" = "Lag phase",
                                                            "cv_area_batch" = "CV cell size",
                                                            "mean_AR_batch" = "Cell shape",
                                                            "mean_area_batch" = "Cell size",
                                                            "var_area_batch" = "Variance cell size")), scales = "free_y") 

response_38_plot

2) Response of the populations that, after adaptation to 38 °C, returned to the control temperature of 20 °C.

# select the populations that experienced 20 °C
response_data_20 <- response_data %>% filter(temperature == 20)

# create variable for time difference since return to 20 °C started
response_data_20 <- response_data_20 %>% mutate(batches_20 = batch - 4)

# create a list with the prior information, using an uninformative prior
# same prior as 38°C models
prior_1_new <- list(R = list(V = 1, nu = 0.002), 
                G = list(G1 = list(V = 1, nu = 0.002)))

# run one separate model for each of the 5 response variables
models_20 <- response_data_20 %>% group_by(data_type) %>% nest() %>%
  mutate(mixed_model = map(data, ~MCMCglmm(fixed = percent_batch1 ~ batches_20, 
                                           random = ~lineage, 
                                           data = ., 
                                           family = "gaussian",
                                           prior = prior_1_new, 
                                           nitt = 2000000,
                                           burnin = 30000, 
                                           thin = 1000,
                                           verbose = F)))

# create new data
new_data_20 <- data.frame(expand.grid(lineage = c("1"), batches_20 = c(0, 1), 
                       percent_batch1 = 0))

# get confidence interval from each model
models_20 <- models_20 %>% mutate(predictions = map(mixed_model, ~predict.MCMCglmm(object = .,
                                                               newdata = new_data_20,
                                                               interval = "confidence")))
models_20 <- models_20 %>% mutate(predictions = map(.x = predictions, ~as_tibble(.)))

# merge new data and model fit for plotting
models_20 <- models_20 %>% mutate(predictions = map(.x = predictions, ~cbind(new_data_20, .))) 

# store models in a separate object
model_20_predictions <- unnest(models_20, predictions)

# plot the data and the model 
response_20_plot <- ggplot(data = response_data_20, aes(x = batch, y = percent_batch1, color = name)) +
  geom_point(aes(color = as.factor(lineage)), alpha = 0.8, size = 2) +
  geom_smooth(data = model_20_predictions,
              aes(x = batches_20 + 4, y = fit, ymin = lwr, ymax = upr), 
              colour = "black", stat = "identity") +
  geom_hline(yintercept  = 0, linetype = "dashed") +
  scale_x_continuous(breaks = c(4, 5)) +
  xlab("Batch") +
  ylab("Change in comparison to control (%)") +
  scale_color_manual(values = c("#d2e0ee", "#90b2d5", "#4d84bc", "#2166ac"),
                     name = "Population", 
                     labels = as_labeller(c( "1" = "Pop. 1",
                                             "2" = "Pop. 2",
                                             "3" = "Pop. 3",
                                             "4" = "Pop. 4"))) +
  theme(text = element_text(size = 15), plot.title = element_text(hjust = 0.5, margin = unit(c(0.5,0.5,0.5,0.5), "lines")),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_blank(),
        legend.position = "bottom") +
  facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
                                                            "lambda" = "Lag phase",
                                                            "cv_area_batch" = "CV cell size",
                                                            "mean_AR_batch" = "Cell shape",
                                                            "mean_area_batch" = "Cell size",
                                                            "var_area_batch" = "Variance cell size")), scales = "free_y") 


response_20_plot

Check the convergence of the MCMC chains.

1) Response to 38 °C models

2) Response to 20 °C models

Fig. 5 Change in population dynamics and morphological traits of T. thermophila populations.

Change in maximum growth rate, lag phase, cell size, variance of cell size, coefficient of variation (CV) of cell size and cell shape are shown, for each population and for each batch separately, using the control cultures at 20 °C in batch 1 as a reference. Maximum growth rate difference is expressed in day-1 and lag phase difference is expressed in days, while all other values are expressed as percent difference. The lines represent the fitted mixed effects models and the shaded areas represent the 95 % credible interval (see methods for details). The colors indicate the temperature in which the population was grown, and the shades represent the population replicate. The dashed lines mark no change in comparison to the control.

# plot both models in the same plot with all the response data

response_plot <- ggplot(data = filter(response_data, data_type != "K"), aes(x = batch, y = percent_batch1, color = name)) +
  geom_point(alpha = 0.8, size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
                     name = "Population and\n temperature", 
                     labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
                                            "2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
                                            "3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
                                            "4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
  geom_smooth(data = model_20_predictions %>% filter(data_type != "K"),
              aes(x = batches_20 + 4, y = fit, ymin = lwr, ymax = upr), alpha = 0.2, colour = "#82a8d0", stat = "identity") +
  geom_smooth(data = model_38_predictions %>% filter(data_type != "K"), 
              aes(x = batches_38 + 2, y = fit, ymin = lwr, ymax = upr), alpha = 0.2, colour = "#d37d88", stat = "identity") +
  geom_hline(yintercept  = 0, linetype = "dashed") +
  xlab("Batch") +
  ylab("Change in comparison to control (%)") +
  theme(text = element_text(size = 15),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.background = element_blank(),
        legend.position = "bottom",
        plot.margin = unit(c(1.5, 1.5, 1.5 ,1.5 ),"cm")) +
  facet_wrap(~data_type, ncol = 3,  labeller = as_labeller(c("mumax" = "Growth rate",
                                                            "lambda" = "Lag phase",
                                                            "cv_area_batch" = "CV cell size",
                                                            "mean_AR_batch" = "Cell shape",
                                                            "mean_area_batch" = "Cell size",
                                                            "var_area_batch" = "Variance cell size")), scales = "free_y") 

response_plot

Table with response estimates and confidence interval for each variable at the final experiment batch.

model_38_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_38 == 3)

## # A tibble: 6 x 7
## # Groups:   data_type [6]
##   data_type              data mixed_model batches_38     fit     lwr    upr
##   <fct>          <list<df[,7> <list>           <int>   <dbl>   <dbl>  <dbl>
## 1 mumax              [20 × 7] <MCMCglmm>           3   4.72    0.468   9.08
## 2 lambda             [20 × 7] <MCMCglmm>           3   0.756  -0.614   2.15
## 3 cv_area_batch      [20 × 7] <MCMCglmm>           3  32.9     5.75   56.6 
## 4 mean_AR_batch      [20 × 7] <MCMCglmm>           3 -16.5   -19.0   -14.1 
## 5 mean_area_bat…     [20 × 7] <MCMCglmm>           3 -30.2   -49.8   -11.2 
## 6 var_area_batch     [20 × 7] <MCMCglmm>           3 -17.4   -39.5     4.79

model_38_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_38 == 0)

## # A tibble: 6 x 7
## # Groups:   data_type [6]
##   data_type                 data mixed_model batches_38    fit    lwr   upr
##   <fct>           <list<df[,7]>> <list>           <int>  <dbl>  <dbl> <dbl>
## 1 mumax                 [20 × 7] <MCMCglmm>           0   1.58  -2.30  5.68
## 2 lambda                [20 × 7] <MCMCglmm>           0   7.72   6.46  9.09
## 3 cv_area_batch         [20 × 7] <MCMCglmm>           0  43.5   19.6  69.2 
## 4 mean_AR_batch         [20 × 7] <MCMCglmm>           0 -10.7  -12.9  -8.38
## 5 mean_area_batch       [20 × 7] <MCMCglmm>           0 -27.3  -47.5  -9.27
## 6 var_area_batch        [20 × 7] <MCMCglmm>           0  13.5   -6.14 36.0

model_20_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_20 == 1)

## # A tibble: 6 x 7
## # Groups:   data_type [6]
##   data_type            data mixed_model batches_20     fit     lwr      upr
##   <fct>         <list<df[,> <list>           <dbl>   <dbl>   <dbl>    <dbl>
## 1 mumax             [8 × 7] <MCMCglmm>           1   0.442  -0.438   1.47  
## 2 lambda            [8 × 7] <MCMCglmm>           1   0.557  -0.607   1.64  
## 3 cv_area_batch     [8 × 7] <MCMCglmm>           1 -10.3   -13.6    -6.62  
## 4 mean_AR_batch     [8 × 7] <MCMCglmm>           1  -5.16  -10.6     0.0956
## 5 mean_area_ba…     [8 × 7] <MCMCglmm>           1 -11.1   -16.6    -5.01  
## 6 var_area_bat…     [8 × 7] <MCMCglmm>           1 -36.9   -43.8   -30.5